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Abstract 


The evolution of unsteady boundary layers on oscillating airfoils is 
Investigated by solving the governing equations by the Characteristic Sox 
scheme. The difficulties associated with coirputing the first profile on a 
given time-line, and the velocity profiles with partial flow reversal are 
solved. A sample calculation has been performed for an external velocity 
distribution typical of those found near the leading edge of thin 
airfoils. The results demonstrate the viability of the calculation 
procedure. 


Introduction 


The effect of unsteady motion of an airfoil on Its stall behavior Is a 
problem of great Importance in many types of fluid motion, for example, 
helicopter rotor blades and jet engine conpressors. A recent detailed 

study has been published by Carr, McAlister and McCroskey^^^ ; this study 
reveals a very complicated phenomenon which depends In a subtle way on a 
large number of parameters. There Is one Important characteristic of the 
data they have compiled which serves as a focal point of the approach 
which we are planning to apply sometime later In our studies; namely that 
at some stage of the cycle, a large vortex Is formed near the surface of 
the airfoil and very shortly afterwards stall occurs. It seems also that 
the stall Is associated with flow reversals In the unsteady boundary layer 
which may spread downstream or upstream de;iend1ng on the leading-edge 
radius of the airfoil. 

The present paper concerns itself with the calculation of boundary -layer 
characteristics of an oscillating airfoil in order to Investigate the 
evolution of unsteady boundary layers on such airfoils. It Is one phase 
of a study which will be extended later in the hooe of throwing light on 
the dynamic stall problem. In our present study we focus our attention on 

*rros£>nt I'll at tlu’ iUTAM Symposium on Unsteady Turbulent Siiear Flows, 
Toulouse, S-8 May l^^Sl. 



the calculation of tlne-depandant boundary laytrs for a givon prtssura 
d1 stributlon. 

There are three difficulties associated with unsteady boundary layers that 
requires careful attention. First of all appropriate Initial values at 
t*0 nust be chosen for the velocity distribution. Strictly speaking they 
can be arbitrary but In that event, the values of 9u/3t at t ■ 0 Is 
non-zero and this Implies an Inviscid acceleration of the fluid In the 
boundary layer and In consequence a velocity of slip begins to grow at the 
wall. This Is smoothed out by an Inner boundary layer Initially of 
thickness (vt)^^^ In which viscous forces are of Importance. Thus a 
double structure develops In the boundary layer which may be treated by a 
generalization of the Keller-box scheme^^^. However, If Interest Is 
centered on the solution at large times, this feature may be rc^duced In 
Importance by requiring that the Initial velociti distribution satisfies 
tne steady-state equation with the Instantaneous external velocity. In 
addition it Is necessary to smooth out the external velocity Ug(x,t) so 
that 3Ug/3t >0 at t - 0 and then standard methods may be used and 
are stable. The use of the smoothing function makes for some loss of 
accuracy at small values of t but the error soon decays to zero once the 
required value of u^ Is specified. In the present problem the choice 
of parameters In the specified external velocity distribution is such that 
the smoothing function Is actually unnecessary but It can easily be 
incorporated Into the scheme. 

The second difficulty arises when u changes sign over part of the 
profile at some x-station where the x-axIs Is parallel to the streamwise 
direction and u Is the corresponding velocity conponent. Normally this 
does not occur and one can Integrate away from the profile In the 
direction of positive u without any difficulty by using a standard 
numerical scheme. However, If the change In sign of u does occur, we 
encounter numerical Instabilities since in the negative u-region we would 
be integrating against the stream. The instability can be avoided by 
Changing the scheme either to the zig-zag box or the characteristics box. 
These new schemes have already been shown to be effective In such 
circumstances when the flow is unsteady^'^ and in three-dimensional 
flows^^"®^. The essence of these schemes is that, to an increasing 
extent, they take into account the fact that small disturbances are 
carried along with the local fluid velocity. 



The third difficulty arises when It Is desired vo compute the first 
velocity profile at the new time line. Given, as we are, the complete 
velocity profile distribution on the previous time line, there Is In 
principle no difficulty In computing values on the next time-line by an 
explicit method, but If we wish to avoid the stability problems associated 
with such a nwthod by using an Implicit method, we are Innedlately faced 
with the problem of generating a starting profile on the new time-line. 

In order to explain the problem further. It Is Instructive to see what 
happens to the stagnation point as a function of time. For this purpose 
let us assume that the external velocity distrlbutloh for an oscillating 
airfoil Is given by the fol^'^wiig function. 


u. 


X ♦ (1 + A sin ut) 


( 1 ) 


(1 ♦ x2)‘s 

where A and 5 ^ denote parameters that need to be specified. This 
equation Is a good approximation to the external velocity distribution 
near the leading edge a class of thin airfoils at variable angles of 
attack and, when A > 0, has recently been used by Cebeci, Stewartson, and 

WilHams^^^ to study leading-edge separation In steady flow. 


Since by def.nitlon 
is given by 


u • 0 at the stagnation point. Its location, x,, 
e 5 


*S “ "^0 ^ 


( 2 ) 


and so the upper and lower surfaces of the airfoil as functions of time 
are defined In particular by 
us take A ■ 1, 

Shown In Figure 1 for one cycle (G ^ t ^ 8). When t ■ 2, the stagnation 


X > x^ and x < x^. For example, let 
w ■ ff/4 and plot Xj/Cq In the (t,x) plane, as 


point x^ is at 


■2 Cq. when t ■ 6 It is at 0, etc. 


If Xj were 


X ■ Xj for all tine and all y. 


fixed we could assume that u * 0 at 
Further profiles at this time-line then follow by use of one of the box 
schemes that have been developed. However x^ i-^ not fixed and it is 
clearly unjustified to assume a priori that u ■ 0 there. Instead 
therefore we use the characteristic box, with an 






Fig. 1. Variation of stagnation 
point with time for one cycle 
according to Eq. (2), with 
u “ ir/4. 


extrapolated normal velocity! to 
complete the profile at x ■ x^, the 
nearest x station to x^ on the 
new time-line. Thus we avoid the need 
to use two x-statlons at the new time- 
line to dimenilonallze the governing 
equations. Once this step Is 
completed, the determination of 
further profiles on the new tine line 
proceeds normally and moreover we can 
Improve our Initial estimate of the 
normal velocity at x ■ x^ by 
Iteration. 


Governing Equations 


The governing boundary-layer equations for an Incompressible laminar or 
turbulent flow past an oscillating airfoil are well known and, with the 
eddy viscosity (C|^) concept, they can be written as 


9u ^ 3v 
dx 3y 


0 


(3) 


u V 

3t 3x 3y 


»Ua ^‘^e 3 t 

3t e ax 3 y 

where t ■ v(3u/3y) - iTV*. In the absence of mass transfer, these 
equations are subject to the boundary conditions given by 


(4) 


y-0, u-v-0; y-^j, u-*- u^(x,t) 


(5) 


The presence of the Reynolds shear stress term, - u' v' , requires a 
closure assumption; In our study we use the algebraic eddy-v1scos1t> 
formulation developed by Cebeci and Smith. For details, see ref. 8. 


To complete the formulation of the problem. Initial conditions must be 
specified in the (t,y) plane at some x “ Xg either on the lower or 
upper su -fac^i of the airfoil (see Fig. 1) as well as initial conditions in 
the (x,y^ plane on both surfaces of the airfoil. In the latter case, if 
we assume that steady-flow conditions prevail at t ■ 0, then the initial 
conditions in the (x,y)-plane can easily be generated for both surfaces by 


solving the governing equations for steady flow, which In this case, are 
given by £q. (3) and by 


u 




du 


e dx 3 y ' ay ' 


(« 5 ) 


where b « v + e_. There Is no problem with the Initial conditions for 

Itl 

Eqs. (1) and (6) since the calculations start at the stagnation point. 


Generation of the Initial conditions In the (t,y) plane at x - x^, 
which Is one ofithe purposes of our study. Is not so easy as was discussed 
In the previous section. The following section describes our solution 
procechjre at t ■ 0 and t > 0. 


Solution Procedure 


We use Keller's Box method to solve the governing equations of the 
previous section. This is a two-point finite-difference method which has 
been used to solve a wide range of parabolic partial-differential 

Q 

equations as discussed by Bradshav, Cebeci and Whitelaw . The solution 
procedure for t ■ 0 (Eqs. (3) and (6)) and t > 0 (Eqs. (3) and (4)) are 
described separately below. 


Solution Procedure for t « 0 

As explained in the Introduction any velocity profile may be chosen at 
t * U to Initiate the computation but It Is convenient to select one 
which obviates the need for double-structured numerical schemes and joins 
smoothly on to the solution for t > 0. We Insist that these profiles 
satisfy the steady state equations with u^ given by (1) and t • 0. 

The details of the procedure for computing these profiles differs slightly 
from previous procedures used In steady two-dimensional flows where we 
have used the definition of the stream function and reduced Eqs. (3) and 
(6) to two first-order ordinary differential equations and to one 
partial -differential equation. Here we consider the solution of Eqs. (3) 
and (6j without the use of the stream function. For this purpose, with 
primes denoting differentiation with respect to y, we let 


-S- 



u' - f 


(7i) 


and write Eqs. (3) and (6) as 


V* 


3u 

3X 


(7b) 


(bf)' 



(7c) 


The finite-difference approximations of Eqs. (7) are also somewhat 
different than those reported In our previous studies dealing with 
two-dimensional flows^®. All quantities except for the normal velocity 
component v, are centered at the center of the box (yj.i/g’ 

*1 1/2^’ *** values of each parameter, say q, at 

the four corners of the Box, that Is, 


1 - 1/2 1 ,^ 1 - 1 / 2 ^ . 1 - 1 /? 


^j-1/2 






'> • * "j'' * “j-i * '■jll* 


However, the centering of the y-velocity component v Is done by writing 
it-as 


1 - 1/2 

j-1/2 


1 


'j-i 


(8b) 


The unknown parameters In Eqs. (8) correspond to Oj and 
that when a solution of the system given by Eqs. (7) Is obtained, f and 
u are computed at (1,j) and v at (1-./2, j). This modified 

centering procedure Is necessary In 
order to avoid oscillations due to the 
use of the continuity equation In the 
form given by Eq. (7b) rather than the 
use of the stream function. The 
centering of Eqs. (7) and the 
subsequent 1 Inearizatlon procedure by 
Newton's method allows the resulting 
linear system to be written In the 


1 - 1/2 
Vj so 



F1g. 2. Net rectang'ia for 
finite-difference approximations. 


form 


( 10 ) 


""j 


6u 


j-1 


h. , 
J-1 


(5f. . dfj.j) 


^^3)j-l 


(9a) 
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(9b) 


(««j - * (v>j («“j * ■ Ci)j 

After llncerlzatlon, the wall boundary conditions 

■ 0, • 0 (10a) 

0 0 

and the edge boundary condition 

Uj ■ (1W>) 

become 

*0, 5v ■ 0, 6u , ■ 0 (11) 

0 0 J 

The resulting linear system consisting of those given in Eos. (9) and by 
those in Eq. (11) can be solved by the block elimination method discussed 
by Cebeci and Bradshaw' 


Solution Procedure ^or t > 0 

We have already pointed out that If au^/at ^ 0 when t ■ 0, a 
double structure scheme should strictly be used to advance the solution 
from t • 0. However the choice of parameters In our study Is such that 
aUg/at is small and the difficulties that arise from using a standard 
method are of a sufficiently minor nature that no further refinement is 
necessary. For larger values of the relevant parameters It is easy to 

Incorporate a smoothing function into u and one can always use the 

( 2 ) • 

general method' . 

Nevertheless there is still the difficulty about obtaining the velocity 
profile on the first x-station at any new time-line. It can be resolved 
with the use of the characteristic box method developed by Cebeci and 
Stewart son?*). Defining the streamline by 

dx 

■ 

u 


dt 

T“ 


- 7 - 


( 12 ) 


and using th« definition of f ' , 
streamline, we write Eq. (4) as 



Fig. 3. Notation and finite- 
difference molecule for the 
Characteristic Box 2. 


and with s denoting the local 



To obtain the solution of the unsteady 
boundary.layer equations given by Eqs. 
(7a,b) and (13) at the first x-statlon 
on either side of the stagnation line, 
let us consider the grid of Fig. 1 and 
direct our attention to the point 
denoted by 1 (see Fig. 3). 


To write the difference approximations of Eq. (13) we define 


ASj • At/cos Oj (14) 

where, with u^ denoting an average velocity, we compute Oj from 

Oj ■ tan"^ Uj (15) 

assuming that at first v at point P is known and Is equal to Its value 
1-1/2 k-1 

at V . This assumption decouples the continuity equation, Eq. 

(7‘ ) from Eqs. (7a) and (13) and reduces the problem to a 
"two-dimensional'* one with f and u being the only unknowns. The 
finite-difference approximations of Eq. (7a) are written In the usual way 
and the finite-difference approximations to Eq. (13) are written by 
centering It at point P. This procedure leads to 


. -1 , -1 ,k -1 ,k. 1 ,k - 


( 16 a) 




. , 1,k n,k-l 

n,k-l u.’ - u.’ 


1 I I - 

^[(Nl+u^) *(Vl*u‘) 3 (■■ 


2 Lr«l + u Jj-i/2 * “ 'j-1/2-' AS. 


(16b) 


- 8 - 


Tht profllts bj* fj «nd «s Mil as at (*i.k-l) art 
obtaintd by Intarpolating the profllts at (1»k-l) and (1-1, k-1). To 
find tht angit m dtfint In Eq. (15) by 

Sinct tht systtm givtn by Eqs. (16) Is lintar, thtrt Is no nttd for 
lintarizatlon and m solvt It subjtct to tht tM boundary conditions, 
namtly. 


“o * “j “ “t 

by using tht block-tllmlnatlon Mthod In which cast tht Mtricts art 
2x2. Wt shall rtftr to this schtnt as Characttrlstic Box 2. 

Once a solution of Eqs. (7a) and (13) has bten obtaintd, wt compute v 
from Eq. (7b) which. In finite-difference form, for the center of the net 
rectangle, point E, can be written as (set Fig. 3) 


V 


i 



1 ,k-l/2 

'Ll 


(19) 


Here denotes the value of at E and u^ Is given by Eq. 

(17). Since the right-hand side of Eq. (19) Is known, m can solve this 

equation for V. and with ■ 0, find V. for 1 < j < J. Wt 
J 0 J — • . p 

then substitute this new value of Vj Into Eq. (16) for v and 

solvt the system again to compute new values of Vj. This procedure Is 
repeated until convergence. 


For convenience we use the same procedure to compute point 2 to the left 
of point 1. Once two points on a given t-llne are cewnputed by this 
procedure, M then use the values of V, at E and E,. compute a 

p J £ 1 

new value v and repeat the solution procedure for Eqs. (7a) and (13). 

and later Eq. (19). After that the stations to the left of point 2 and 

the stat’ons to the right of point 1 are computed by using the Regular Box 

scheme If there Is no flow reversal across the layer and by the Character. 

i*t1c Box scheme If there is flow reversal. The "new" Characteristic Box 

scheme is now slightly different than the Characteristic Box 2 so we shall 

refer to it as Characteristic Box 3. 


- 9 - 


To doscribt tho Charocttrlstic Box 3 sch«M which solvot Eqs. (7a. b) and 
(13) without dtcoupling th« continuity t^atlon fnoai Eqs. (7a) and (13). 
wt considtr tht skttch shown In Flqurt 4. 


' 1 - 1/2 


'1-5/2 


2 l 


1 


V • 

''1-3/2 


\ 








1 


n c 




1-2 1-1 ' 1 1-»*1 1+2 
Fig. 4. Notation for Characttrlstic Box 3. 

using tha Zig-Zag Box schema discussed In detail In Ref. 5. we write Eq. 
(7b) In the following finite-difference form 


^*-1 "■^'' j - 1/2 * 'h * ^ j - 1/2 


where 


e 


' ^ 1 ^ 

(* 1+1 ■ * 1 .l} ■ *<_l) 

X^ - X 

h • X 


'1 


‘ 1-1 


‘ 1+1 ’ * 1-1 ^ 


‘1 ■ '‘ 1 -r 

i 


* i+J 


'“j-1/2 “j-1/2 ' 


Since 


.y , ^^- 1/2 • ^-3/2| . 

1-1/2 ''1-3/2 (* 1 . 3/2 - * 1 . 5 / 2 } ''1-5/2^* 


3/2 '‘1-5/2' 

The relation between v^.j /2 **• *»'‘^Tten as 


* V ^ 1 - l /2 * ®3 

where 62 83 *'*® 


( 20 ) 


( 21 ) 


( 22 ) 


®2 “ '1*7: 


^*< " '‘ 1 - 1 / 2 ^ 

3/2 ■ *1-1/2^ 


*3 * 'Vl-3/2 


Introducing Eq. (22) Into Eq. (20) and rearranging, we get 
.-1 


where 


h]-i (1 - *^«j-l/2 ^ 84 

64 ■ -81 ♦ «“jli/2'' * ^jii ^83 ^ " 83^ 


(23) 


(24) 


(25) 


- 10 - 



As In Ch«rtct«r1st1c Box 2, Mt ctnttr Eq. (13) «t thq midpoint of (1,k) 
end (l.n) to got tht f1nltt-d1ff«rcnct tguetlons givtn by Eq. (16b) with 
boing obtaintd by llnoir Inttrpolatlon of • 

which It 


• Vj 




(26) 


Equations (16a), (24) and (16b) art than lintarittd by Ntwton's mathod, 
and again art solvad by tha block allmlnatlon mathod. 


Whan thara Is no flow ravarsal across tha layar, wa usa tha Ragular Box 
schoM dascribad In datall In Raf> 5. 


Rasu 1 tt and Discussion 


To data calculations have baan carriad out In only ona tast-casa, namaly 
whan Cq * A - 1, w ■ x/4, and for a timitad ranga of x 

( |x| < 0.3). With tha usa of tha various procedures described In this 
paper the calculations were quite straightforward and tha formal validity 
and efficacy of the numerical schemas were established. Tha results are 
suimnarized In Figs. 5-7. In Fig. 5 wa display tha variation of wall-shear 
with time at different x-statlons and In Fig. 6 tha variation with x at 
different times. These graphs are entirely in line with expectations and 
wa note that tha flow reversal at tha wa1 : Is quite smooth. A similar 
remark applies to tha velocity profiles on either side of tha stagnation 
line displayed In Fig. 7. 

Tha next phase In our studies Is to extend tha computations to larger 

values of A.(^ and smaller values of w to as to more closely 

approach tha conditions of dynamic stall. It Is of Interest to comment on 

the ^luid mechanical problems that may than arise. First, if Is 

increased beyond 1.155, tha steady-state solution at t > 0 separates on 

the upper side of tha airfoil and tha calculation terminates. This is not 

a serious drawback unless tha unsteady boundary develops a singularity 

because the smoothing function mentioned earlier may be adapted to ensure 

that £ Is initially less than 1.155 and rises to a value g*‘eater than 
'0 

that after a finite time. The unsteady boundary layer then includes 
regions of reversed flow which may well become extensive if A is also 


- 11 - 


allowed to increase to mimic more closely the conditions of dynamic 

stall. Even If the boundary layer remains smooth, the displacement 

thickness may then become much thicker and have a significant modifying 

effect on the external flow. It would be useful then to consider an 

Interactive problem wMch the external stream depends in part on the 

displacement thickness thus generalizing the studies reported for steady 

flow^^^. Consideration then has to be given to the variation of 

circulation with time which may lead to a more complicated expression for 

the dependence of u. on the displacement thickness than was used* in 

e 

Ref. 7 but the computation should not be any more complicated as a 

result. Finally, in order to mimic the dynamic stall problem most 

effect! vely^^’^^^ , w should be reduced to very small values (as 

typical of dynamic stall problem). So long as u > 0, the difficulties 

reported in Ref. 3 at separation in uninteracted flows and in 

post-separation flows otherwise, should not be present. On the other 

(121 

hand, van Donsnelen and Shen' ' have 
provided quite strong evidence that a 
singularity can occur in an unsteady 
boundary layer for which the external 
velocity is steady. This phenomenon 
is still somewhat controversial 
but there seems no doubt that the 
boundary layer will exhibit dramatic 
properties for small enough values of 
- 0.2 - C .1 0 0 ^ 0.2 0.2 u and it is possible that these may 

Ftg. 5. of Mil shear 

parameters with time at differ- stall, 
ent x-stations. 




Fig. 6. Variation of skin-friction coefficient with x at different 
t-intervals. Note t ■ 8 corresponds to one cycle. 
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t * 0 



t » 4 


LOI.'ER 


UP^FR 





X = -n.l4 



LOIJFR 


UPPER 



Fig. 7. Velocltr profiles in the immediate neighborhood of the stagnation 
lines at oifferent times, u ■ 0 on the dashed line at the specified time 
and serves to "define" upper*or lower surfaces. 
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